####################################################
####################################################
####################################################

dists <- readRDS("temp/shooting_demos.rds") |>
  mutate(date = as.Date(date),
         id = id4, ## id (and distance, in next line) of nearest 4-victim shooting
         dist = dist4) |>
  filter(!is.na(id),
         id %in% (readRDS("temp/geocoded_shootings_from_gva.rds") |>
                    filter(home + dv == 0))$id ## drop if in the home or domestic violence
  )

### loop over thresholds for RDITS
run_rdd <- function(threshold){
  library(data.table)
  library(tidyverse)
  library(rdrobust)
  print(threshold)
  ##keep the observations within the threshold closest to election day
  ## using data.table notation to speed this up
  set_pre <- dists[dists$dist <= threshold & dists$pre,
                   .SD[date %in% max(date)], by=list(year, GEOID)]
  set_pre <- set_pre[, .SD[1], by = .(year, GEOID)] %>% 
    mutate(treated = T)
  
  ##keep the observations within the threshold closest to election day
  set_post <- dists[dists$dist <= threshold & !dists$pre,
                    .SD[date %in% min(date)], by=list(year, GEOID)]
  set_post <- set_post[, .SD[1], by = .(year, GEOID)] %>% 
    mutate(treated = F)
  
  full_treat <- bind_rows(set_pre, set_post)
  
  full_treat <- full_treat[, if (.N == 1) .SD, by = .(year, GEOID)]
  
  full_treat <- full_treat |> 
    select(GEOID, id, date, dist, year, turnout,
           median_income, nh_white, nh_black, median_age, pop_dens, latino, asian,
           some_college, turnout_pre, treated, population, av_temp, cases, share_group,
           in_person_rate, contact, homicide) %>% 
    mutate(d2 = ifelse(year == "2020", as.integer(date - as.Date("2020-11-03")),
                       as.integer(date - as.Date("2016-11-08"))),
           t16 = year == "2016")
  
  if(threshold == 0.5){
    saveRDS(full_treat, "temp/primary_half_only_one.rds")
  }
  ########################################
  
  
  l <- rdrobust(y = full_treat$turnout, x = full_treat$d2, p = 1, c = 0, cluster = full_treat$id,
                covs = select(full_treat,
                              latino, nh_white, asian,
                              nh_black, median_income, median_age,
                              pop_dens, turnout_pre, t16, some_college))
  

  
  f <- tibble(coef = l$coef,
              se = l$se, 
              pv = l$pv,
              p = threshold,
              n_pre = l$N_h[1],
              n_post = l$N_h[2],
              n = l$N_h[1] + l$N_h[2],
              bw = l[["bws"]][1],
              u = l[["ci"]][,2],
              l = l[["ci"]][,1])
}

cl <- makeCluster(8)  
registerDoParallel(cl)

clusterExport(cl, list("run_rdd", "dists"))
out <- rbindlist(parLapply(cl, c(0.25, 0.5, 0.75, seq(1, 10, 1)), fun = run_rdd))

out$estimate <- rep(c('Traditional','Bias-Adjusted','Robust'),nrow(out)/3)

out <- mutate_at(out, vars(coef, l, u), ~. * -1)

out$estimate <- factor(out$estimate, levels = c('Traditional','Bias-Adjusted','Robust'))

saveRDS(out, "temp/only_one_out_data.rds")

out <- readRDS("temp/only_one_out_data.rds")

different_dists <- ggplot(filter(out, estimate == "Robust", p <= 10),
                          aes(x = p, y = coef, ymin = l, ymax = u)) +
  geom_point() +
  geom_errorbar(width = 0) + 
  theme_bc(base_family = "LM Roman 10", base_size = 10,
           plot.margin = margin()) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_x_continuous(breaks = seq(0, 10, 2)) +
  scale_y_continuous(limits = c(-0.1, 0.2)) +
  labs(y = "Local average treatment effect", x = "Radius around shooting (miles)")
different_dists

saveRDS(different_dists, "temp/different_dists_solo_obs.rds")

for(i in outpaths){
  ggsave(paste0(i, "solo_obs_results.png"), height = 4, width = 3, units = "in")
}
